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Abstract 

We report Monte Carlo results for the two-dimensional hard disk system in the 
transition region. Simulations were performed in the NVT ensemble with up to 
1024 2 disks. The scaling behaviour of the positional and bond-orientational order 
parameter as well as the positional correlation length prove the existence of a hex- 
atic phase as predicted by the Kosterlitz-Thouless-Halperin-Nelson- Young theory. 
The analysis of the pressure shows that this phase is outside a possible first-order 
transition. 
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The nature of the two-dimensional melting transition has been an unsolved 
problem for many years [1,2]. The Kosterlitz-Thouless-Halperin- Nelson- Young 
(KTHNY) theory [3,4,5,6] predicts two continuous transitions. The first tran- 
sition occurs when the solid (quasi-long-range positional order, long-range ori- 
entational order) undergoes a dislocation unbinding transition to the hexatic 
phase (short-range positional order, quasi-long-range orientational order). The 
second transition is the disclination unbinding transition which transforms this 
hexatic phase into an isotropic phase (short-range positional and orientational 
order). An alternative scenario has been proposed by Chui [7,8]. He presented 
a theory via spontaneous generation of grain boundaries, i.e. collective excita- 
tions of dislocations. He found that grain boundaries may be generated before 
the dislocations unbind if the core energy of dislocations is sufficiently small, 
and predicted a first-order transition. This is characterized by a coexistence 
region of the solid and isotropic phase, while no hexatic phase exists. Another 
proposal was given by Glaser and Clark [2]. They considered a detailed the- 
ory where the transition is handled as a condensation of localized, thermally 
generated geometrical defects and found also a first-order transition. Calcu- 
lations based on the density-functional approach were done by Ryzhov and 



Preprint submitted to Elsevier Science 



2 February 2008 



Tareyeva [9]. They derived that the hexatic phase cannot exist in the hard 
disk system. 

Even for the simple hard disk system no consensus about the existence of 
a hexatic phase has been established thus far. The melting transition of the 
hard disk system was first seen in a computer simulation by Alder and Wain- 
wright [10]. They used a system of 870 disks and molecular dynamics meth- 
ods (NVE ensemble) and found that this system undergoes a first-order phase 
transition. But the results of such small systems are affected by large finite-size 
effects. Simulations performed in the last years used Monte Carlo (MC) tech- 
niques either with constant volume (NVT ensemble) [11,12,13,14,15,16,17,18,19,20,21] 
or constant pressure (NpT ensemble) [22,23,24]. Zollweg, Chester and Le- 
ung [11] made detailed investigations of large systems up to 16384 particles, 
but draw no conclusives about the order of the phase transition. The anal- 
ysis of Zollweg and Chester [12] for the pressure gave an upper limit for a 
first-order phase transition, but is compatible with all other scenarios. Lee 
and Strandburg [22] used isobaric MC simulations and found a double-peaked 
structure in the volume distribution. Lee-Kosterlitz scaling led them to con- 
clude that the phase transition is of first-order. However, the data are not in 
the scaling region, since their largest system contained only 400 particles. MC 
investigations of the bond-orientational order parameter via finite-size scaling 
(FSS) with the block analysis technique of 16384 particle systems were done 
by Weber, Marx and Binder [13,14]. They also favoured a first-order phase 
transition. In contrast to this, Fernandez, Alonso and Stankiewicz [23,24] 1 
predicted a one-stage continuous melting transition, i.e. a scenario with a sin- 
gle continuous transition and consequently without a hexatic phase. Their 
conclusions were based on the examination of the bond-orientational order 
parameter of different systems up to 15876 particles and hard-crystalline wall 
boundary conditions. Mitus, Weber and Marx [15] studied the local structure 
of a system with 4096 hard disks. From the linear behaviour of a local order 
parameter they derived bounds for a possible coexistence region. Sengupta, 
Nielaba and Binder [16] simulated a dislocation free triangular solid of hard 
disks using a constrained Monte Carlo algorithm and showed that a KTHNY 
transition preempts a first-order transition. Combining renormalisation groups 
ideas with MC input they derived also an estimate of p m = 0.914(2) for a pos- 
sible hexatic-to-crystal transition [17]. Finally, in [18,19] the present author 
determined the disclination binding transition density p; = 0.899(1) from mea- 
surements of the bond-orientational correlation length and susceptibility in the 
isotropic phase as well as the critical exponent r] 6 = 0.25(4) from finite-size 
scaling of systems with up to 16384 particles. The results are in agreement 
with the KTHNY theory, while a first-order phase transition with small cor- 
relation length and a one-stage continuous transition can be ruled out. By 
studying short-time behaviour and FSS of the positional order [20,21] we also 
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estimated the dislocation binding density with p m 0.933 as well as the 
critical exponent rj xs 0.200. 

In this letter, we present results obtained through MC simulations in the 
NVT ensemble to answer the question of the existence of a hexatic phase 
and therefore the kind of the phase transition. Although we have already 
shown that the disclination binding density p x = 0.899(1) is lower than the 
dislocation binding density p m ~ 0.933, no direct observation of the hexatic 
phase so far exists. Therefore, we examine the positional order as well as the 
orientational order within this region at p = 0.914 and p = 0.918. Simulations 
in the NVT ensemble are computationally less expensive than in the NpT 
ensemble. The density region is outside a possible first order transition thus 
there is no necessity to perform constant pressure simulations. 

We consider systems of N = 32 2 up to 1024 2 hard disks in a two-dimensional 
rectangular box with ratio v^3 : 2. Additionally, we studied three systems 
in a square box to determine the influence of boundary conditions. The disk 
diameter is set equal to one. For the simulations we used an improved updating 
scheme [27], in which the conventional Metropolis step of a single particle is 
replaced by a collective (non-local) step of a chain of particles. 

Due to the improved updating scheme, which reduces autocorrelation times, 
we were able to perform simulations of large systems within the transition 
region. Although, the scaling behaviour is probably not changed compared 
to the Metropolis algorithm, the chain algorithm significantly speeds up the 
simulations. Especially in the hexatic phase with large positional order length, 
this algorithm is advantageous 2 . 

We measured the kth moment of the global bond-orientational order param- 
eter 

k 

(1) 



where the sum on j is over the iVj neighbours of this particle and 6ij is the 
angle between the particles % and j and an arbitrary but fixed reference axis. 
Two particles are defined as neighbours if the distance is less than 1.4a, where 

a = \jlj\pip is the average lattice spacing of a perfect crystal. This definition 
is computationally less expensive than precise determination by the Voronoi 
construction. While the values of ip 6 shown in table 1 depend on the definition 
of neighbours, the conclusions are independent of the method used. 



2 The chains or particles moved in a single step are much longer than those observed 
for p = 0.898 [27]. 
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The fcth moment of the positional order parameter is defined as 

k 

(2) 



where G is a reciprocal lattice vector and denotes the position of particle i. 
The magnitude of G is given by 2n/a, while the orientation was determined 
from the global bond-orientation. With ip 6 ~ exp(6i#) we defined the angle 
9 (— 7r/6 < 9 < 7r/6) and therefore the orientation of the crystal and G. The 
fourth-order cumulant for the positional order is given by 




ib k = 

Ypos 



1 N 

^exp(iG • Ti) 



N 



i=i 



The bond-orientational correlation length £ pos was extracted from the 'zero- 
momentum' correlation function of -?/>pos( x ) 

9 P os(x) = (^-J2^pos(x,y) ) j ^£^pos(0,y')j > ( 4 ) 



by fitting the data with a single cosh, where denotes the number of particles 
in a stripe between x + Ax/2 and x — Ax/2 [19]. The pressure is calculated 
from the pair correlation function g(r) 

■P l + -ps(l) , (5) 



NkT 2 ' V 2 



where A is the close-packed area of the system, i.e. A = N\/3/2. 

All simulations of systems in the rectangular box started from the ordered 
state while we used a closed-packed state as initial condition in case of the 
square box 3 . Careful attention has been paid to the equilibration of the sys- 
tems. For that purpose we estimated autocorrelation times from binning 4 . 

3 Disordered initial states for the hard disk model at high densities are more difficult 
to simulate and equilibration is more time consuming. 

4 We built blocks of subsequent configurations, called bins, and averaged the quan- 
tities first in the bin. The obtained bin averages themselves can be considered as the 
results of single measurements. If the bins are large enough, then the average values 
in different bins are practically uncorrelated and the obtained statistical errors are 
correct. All given statistical errors are obtained in this way, i.e. taking correlations 
into account. This is also a way to estimate the autocorrelation time. 
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Table 1 

The density of three and eight coordinated particles, the pressure, the second mo- 
ment of the global bond-orientational and positional order parameter, the positional 
fourth-order cumulant and the positional correlation length for different densities, 
boundary conditions and system sizes. 
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Table 2 

Estimated integrated autocorrelation time measured in units of thousand sweeps 
of the chain Metropolis algorithm. We used different (optimized) step sizes for each 
density. 
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Additionally, two systems were studied in detail, i.e. we monitored all quanti- 
ties and measured the autocorrelation times. We used the defect density as an 
estimate for the equilibration of dislocations. The results are given in table 2. 
The autocorrelation times are much smaller than the number of warm-up 
sweeps, which were at least one million. The number of measurement sweeps 
ranged from 2 to 20 millions. 

For a hexatic phase — as predicted by the KTHNY theory — FSS implies 
ipQ 2 ~ L~ V6 for the orientational order parameter, while the positional order 
parameter should scale with il) po 2 ~ L~ 2 for large enough system sizes (L ^> 
£ pos )- The cumulant U pos should decrease with increasing system size L for 

P < Pin- 



We measure ?/> 6 2 , VVs 2 > U pos , £ pos and pA /NKT as well as the density of 
three (n 3 ) and eight coordinated particles (n 8 ) 5 . We use systems of N = 
32 2 ,64 2 , 128 2 ,256 2 ,512 2 and 1024 2 particles. All results are given in table 1. 
Statistical errors have been calculated by binning. We also added systematic 
errors coming from the interpolation to g(l) for the pressure as well as errors 
coming from the interval used for fitting in case of the correlation length. 
Simulations of the hard disk system are not only affected by usual finite-size 
effects, but also by systematic errors coming from the boundary conditions. 
Even for a rectangular box of ratio v^3 : 2 and (quasi-) long-range order crystal 
tilting occurs. This leads to large autocorrelation times as well as additional, 
complicated finite-size effects compared to simple lattice models. 

The pressure and positional correlation length at p = 0.918 are nearly inde- 
pendent for system sizes with N > 64. The results of n 3 and n 8 for iV = 32 2 
show the suppression of defects for systems that are too small. Leaving out 
the data for iV = 32 2 the results at p = 0.918 seem to be mostly consistent. 
The results for ipQ and ip pos are shown in fig. 1. Obviously, the orientational 
order is (quasi-) long-range, while the positional order is short-ranged 6 which 
is confirmed by the decrease of U pos . Deviations are caused by statistical er- 
rors coming from large autocorrelation times and systematic errors due to 
boundary conditions. For example, the increase of ipQ 2 for lower system sizes 
(N = 32 2 ) and rectangular boundary conditions is caused by stabilization 
effects. The (quasi-) long-range of the orientational order is consistent with 
previous measurements [19] where p x = 0.899(1) was determined. The mea- 
surements at p = 0.914 confirm the short-range positional order in this range. 

Although, the results should be independent from the initial conditions, our 
data can differ due to the different boundary conditions used. However, the 
data should converge for increasing system sizes. Also, for larger systems the 
correlation length, which is direction dependent for smaller systems, should 
be isotropic. 

A visualisation of the orientational order for a typical configuration at p = 
0.918 with N = 1024 2 particles in the rectangular box is shown in fig. 2a and 
b. We divided the system into 380 2 subsystems with approximately 7 particles 
per subsystem and calculated ipQ k averaged over the particles inside. The left 
picture visualises the local orientation. The initial state of a perfect ordered 
crystal (9 = 0) corresponds to white areas. Areas with \9\ > are grey, where 

5 These values are obtained from the definition for neighbours given above. 

6 We carefully tried to equilibrate the systems. The measurements of the correlation 
time shows that the only critical point is N = 1024 2 . However, even a too short 
equilibration time for the largest system wouldn't affect our conclusions. The reason 
is that due to the initial ordered state for the rectangular boxes, a non-equilibrated 
system would lead to too high values for the order parameter, i.e. the decrease of 
V> P os 2 for increasing system sizes would be even larger. 
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Fig. 1. Double logarithmic plots of the second moment of the bond-orientational 
(left figure) and the positional order (right figure) parameter at p = 0.918 as a 
function of y/N, which is proportional to the system size L. Full symbols are for 
the rectangular box while open symbols denote square boxes. The lines are guides 
for the eye. 
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Fig. 2. Visualisations of the local bond-orientational order parameter ipQ of a typ- 
ical configuration with 1024 2 particles in a rectangular box. In the left picture the 
amount of black represents the local orientation. White areas correspond to the ini- 
tal orientation of the perfect ordered system (9 = 0). The larger the absolute value 
of 0, the darker the areas. The right picture shows the degree of order, i.e. the local 
value of < ip6 2 < 1 is visualised. White areas correspond to perfect ordered areas. 

the amount of black is proportional \9\. The case \9\ = tt/6 corresponds to 
black regions. The right picture shows the local order, i.e. the amount of black 
is chosen proportional to 1 — tpe 2 . Thus, the initial state of a perfect ordered 
crystal corresponds to a white area. The two pictures show that areas with 
orientations significant different from 9 = are normally small, i.e. have less 
than 100 particles. 

Our simulations of the hard disk model in the NVT ensemble at p = 0.914 
and p = 0.918 prove the short-range of the positional order and therefore 
the existence of a hexatic phase as predicted by the KTHNY theory. The 
positional correlation length is about 15 at p = 0.914 and 40 at p = 0.918, 
respectively. The orientational order is quasi-long-range. The scaling of 
yields i] 6 0.015 at p = 0.914 and i] 6 = 0.005(3) at p = 0.918, respectively. 
However, we cannot rule out long-range orientational order. The observed 
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phase is not within a possible first-order phase transition since the pressure 
8.04(1) at p = 0.918 is higher than that of such a transition 7.95(1) [19]. 
Therefore, a one-stage continuous transition [23,24] as well as a first-order 
phase transition from the isotropic to the solid phase can be ruled out. Taking 
the results of previous measurements [18,19] a KTHNY-like phase transition 
is most likely. However, a first-order phase transition from the isotropic to the 
hexatic phase with very large orientational order correlation lengths cannot be 
ruled out. Detailed investigations of the pressure around p- x would be necessary 
to make a decision. Also, the exact value of p m and therefore of rj as well as 
the behaviour of rj e — which should approach zero according to the KTHNY 
theory — have to be examined in order to reach a final conclusion as to the 
question of the kind of transition. 
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